// STL
#include <vector>                      // std::vector<>

// Boost
#include <boost/math/distributions/non_central_chi_squared.hpp> // boost::math::non_central_chi_squared, ::cdf(), ::quantile()

// stats++
#include "statsxx/postprocess/ROC.hpp" // ROC_pt, binormal_to_ROC()


//
// DESC: Calculate (explicitly) the binormal-LR ROC curve.
//
// INPUT:
//
//     double                             d_a   : Metz and Pan binormal-LR parameter
//     double                             c     : Metz and Pan binormal-LR parameter
//
//     unsigned int                       n     : number of points to calculate binormal ROC curve over (default of 100).
//                                                NOTE: The actual number of points stored in the ROC curve will actually be two less, because the end points are discarded.
//
// OUTPUT:
//
//     std::vector<ROC_pt>                ROC   : ROC curve
//
// -----
//
// NOTE: Scores are not calculated for the (returned) ROC curve.
//
// -----
//
// NOTE: See notes in fit_binormal_LR() for relevant formulas and explanations.
//
// -----
//
// NOTE: See:
//
//    C. E. Metz and X. C. Pan, ``"Proper" binormal ROC curves: theory and maximum-likelihood estimation,'' J. Math. Psych. 43, 1-33 (1999).
//
//    S. L. Hillis, ``Equivalence of binormal likelihood-ratio and bi-chi-squared ROC curve models,'' Stat. Med. 35, 2031--2057 (2016).
//
inline std::vector<ROC_pt> binormal_LR_to_ROC(
                                              const double        d_a, // Metz and Pan binormal-LR parameter
                                              const double        c,   // Metz and Pan binormal-LR parameter
                                              // -----
                                              const unsigned int  n
                                              )
{
    std::vector<ROC_pt> ROC;

    //=========================================================
    // (PRE-)CHECK
    //=========================================================

    // NOTE: If c == 0, a proper binormal ROC curve results.
    //
    // NOTE: ... So just call that:
    if( c == 0. )
    {
        return binormal_to_ROC(
                               d_a, // a == d_a
                               1.,  // b == 1.
                               // -----
                               n
                               );
    }

    //=========================================================
    // INITIALIZE
    //=========================================================

    // CONVERT Metz AND Pan PARAMETERS TO Hillis
    //
    // NOTE: See Section 5 in Hillis.

    double lambda = (1. - c)*(1. - c)/((c + 1.)*(c + 1.));

    double theta  = (1./(16.*c*c))*d_a*d_a*(c*c + 1.)*(c + 1.)*(c + 1.);

    double lambda_theta = lambda*theta;

    // SETUP CHI^2 DISTRIBUTIONS
    //
    // NOTE: See Section 4.1.1 in Hillis.

    boost::math::non_central_chi_squared chi2_dist_lambda_theta(1, lambda_theta);
    boost::math::non_central_chi_squared chi2_dist_theta(1, theta);

    //=========================================================
    // CALCULATE BINORMAL-LR CURVE
    //=========================================================

    // NOTE: See the NOTEs in ROC.hpp, related to / justifying plotting uniform in FPR.

    //---------------------------------------------------------
    // CALCULATE FPR AND TPR
    //---------------------------------------------------------
    std::vector<double> FPR;
    std::vector<double> TPR;

    for( unsigned int i = 1; i < (n-1); ++i )
    {
        double _FPR = static_cast<double>(i)/static_cast<double>(n-1);

        // NOTE: See Section 4.1 in Hillis.
        double _TPR;

        if( lambda > 1. )
        {
            // Eq. (13)
            _TPR = 1. - boost::math::cdf(chi2_dist_lambda_theta, ((1./lambda)*boost::math::quantile(chi2_dist_theta, (1. - _FPR))));
        }
        // NOTE: lambda == 1 does NOT occur; for c == 0, binormal_to_ROC() is called. See above.
        else // ( lambda < 1. )
        {
            // below Eq. (14)
            _TPR = boost::math::cdf(chi2_dist_lambda_theta, ((1./lambda)*boost::math::quantile(chi2_dist_theta, _FPR)));
        }

        FPR.push_back(_FPR);
        TPR.push_back(_TPR);
    }

    //=========================================================
    // CONSTRUCT ROC CURVE
    //=========================================================

    for( auto i = 0; i < FPR.size(); ++i )
    {
        ROC_pt pt;
        pt.FPR = FPR[i];
        pt.TPR = TPR[i];

        ROC.push_back(pt);
    }

    // ADD IN THE ASYMPTOTIC POINTS
    ROC_pt pt;

    pt.FPR = 0.;
    pt.TPR = 0.;

    ROC.insert(ROC.begin(), pt);

    pt.FPR = 1.;
    pt.TPR = 1.;

    ROC.push_back(pt);

    //=========================================================
    // RETURN
    //=========================================================

    return ROC;
}
